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Abstract. We study a very simple model of a short-range attraction and an outer 
shell repulsion as a test system for demixing phase transition and density anomaly. The 
phase-diagram is obtained by applying mean field analysis and Monte Carlo simulations 
to a two dimensional lattice gas with nearest-neighbors attraction and next-nearest- 
neighbors repulsion (the outer shell). Two liquid phases and density anomaly are found. 
The coexistence line between these two liquid phases meets a critical line between the 
fluid and the low density liquid at a tricritical point. The line of maximum density 
emerges in the vicinity of the tricritical point, close to the demixing transition. 
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1. Introduction 

Water is anomalous substance in many respects. 

Most liquids contract upon cooling. This is not the case of water, a liquid where 
the specific volume at ambient pressure starts to increase when cooled below T = 4°C 
pp. Besides, in a certain range of pressures, also exhibits an anomalous increase of 
compressibility and specific heat upon cooling [21- It is less well known that water 
diffuses faster under pressure at very high densities and at very low temperatures [H]. 
Besides, the viscosity of water decreases upon increasing pressure 0-11] • 

It was proposed a few years ago that these anomalies are related to a second critical 
point between two liquid phases, a low density liquid (LDL) and a high density liquid 
(HDL) p. This critical point was discovered by computer simulations. This work 
suggests that this critical point is located at the supercooled region beyond the fine 
of homogeneous nucleation and thus cannot be experimentally measured. Even if this 
limitation, this hypothesis has been supported by indirect experimental results [TjpTl]. 

In spite of the limit of 2357^" below which water cannot be found in the liquid 
phase without crystallization, two amorphous phases were observed at much lower 
temperatures [TT]. There is evidence, although yet under test, that these two amorphous 
phases are related to two liquid phases in fluid water [T2] [IHl • 

Water is not an isolated case. There are also other examples of tetrahedrally bonded 
molecular liquids such as phosphorus [Hj [T3] and amorphous silica [1^] that also are good 
candidates for having two liquid phases. Moreover, other materials such as liquid metals 
[T7] and graphite JH] also exhibit thermodynamic anomalies. Unfortunately a coherent 
and general interpretation of the low density liquid and high density liquid phases is 
still missing. 

What type of potential would be appropriated for describing the tetrahedrally 
bonded molecular liquids? Directional interactions are certainly an important ingredient 
in obtaining a quantitative predictions for network-forming liquids like water. However, 
the models that are obtained from that approach are too complicated, being impossible 
to go beyond mean fleld analysis. Isotropic models became the simplest framework to 
understand the physics of the hquid-liquid phase transition and liquid state anomalies. 

Recently it has been shown that the presence of two liquid phases can be associated 
with a potential with an attractive part and two characteristic short-range repulsive 
distances. The smallest of these two distances is associated with the hard core of the 
molecule, while the largest one with the soft core pC^^ Acknowledging that core softed 
(CS) potentials may engender a demixing transition between two liquids of different 
densities, a number of CS potentials were proposed to model the anisotropic systems 
described above. The flrst suggestion was made many years ago by Stell and coworkers in 
order to explain the isostructural solid-solid transition ending in a critical point jSOj- [211 • 
Debenedetti et al. using general thermodynamic arguments, confirmed that a CS can 
lead to a coefficient of thermal expansion negative and consequently to density anomaly 
[22] • This together with the increase of the thermal compressibility has been used as 
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indications of the presence of two liquid phases j2Sl|211 which may be hidden beyond 
the homogeneous nucleation. The difficulty with these approaches is that continuous 
potentials usually lead to crystallization at the region where the critical point would be 
expected. The analysis of the presence of both the two liquid phases and the critical 
point becomes indirect. 

Nevertheless the apparent success of the CS potentials, it is not clear that the 
presence of two liquid phases and anomalies are associated only with them [22] or if 
they result from the presence of two competing scales like the ones present in the CS 
potentials. In this work we analyze another type of model system where two competing 
scales are also present. We study a potential with a hard core, a short-range attractive 
part and a repulsive shell. While the attraction accounts both for the van der Waals 
and hydrogen bonding interactions, the outer shell repulsion is related to the interstitial 
molecule that break the tetrahedral structure like the interstitial oxygens in water. In 
order to circumvent crystallization, our model system is a two dimensional lattice gas 
with first-neighbors attraction and second neighbors repulsion. The system is in contact 
with a reservoir of particles. 

We show that this very simple system exhibits both density anomaly and two 
liquid phases. However, instead of having a critical point ending the coexistence line 
between the two liquid phases as one usually would expect, it has a tri critical point. 
The connection between the presence of criticality and the density anomaly shown. 

The reminder of the paper goes as follows. In sec© the model is presented, the 
mean field analysis is shown on sec.(jni), results from simulations are discussed in sec.(0]) 
and our findings are summarized in sec.(0). 

2. The model 

Our system is represented by a square lattice with sites. Associated to each site there 
is an occupational variable, cTj. If the site is occupied by molecule, CTj = 1, otherwise 
CTj = 0. Each site interacts with its nearest neighbors with attractive interactions and 
with its next-nearest neighbors with repulsive interactions ( see Fig^. Therefore the 
Hamiltonian of this system is given by 

Ti- = -Vi ai(Tj - Va Yl ^^^^ - ^^Y^i^ (1) 

<ij> «ik» i 

where the first sum is over the four nearest-neighbors, with energy gain Vi > 0, and 
the second sum, over the four next-nearest-neighbors, has energy interaction V2 < 
(see Fig|21). The last term in eq.(0) refers to chemical potential contribution /i over all 
molecules. Here we will consider periodic boundary conditions. 

In order to help visualization of the possible phases, the lattice is divided into four 
sub-lattices: 1,2,3,4 (appropriated for the description of an arbitrary superstructure 
with twice the lattice spacing of the original lattice). The corners of a simple square are 
labeled counter-clockwise to indicate the sub-lattices (see example in Fig 01) 
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Figure 1. Intermolecular interaction 
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Figure 2. The neighbor's interactions 

The sub-lattice density is given by 

P/3 = ^E^.> /? = !,. ..,4. (2) 

The ground state is defined by the lower grand potential free energy at T = 0, and 
therefore the minimum value of Hamiltonian as given by eq.(^. Due the symmetry of 
the system, in the case Vi > and V2 < 0, there are two possibilities illustrated below. 

2.0.1. If Vi < — 2V2 For > — 3Vi — 4V2 the system stands in the the dense hquid 
phase (DL), where all sites are occupied. As the chemical potential is decreased, at 
/i = — 3Vi — 4V2 there is a phase transition from the dense liquid phase to the structured 
dilute liquid phase (SDL), illustrated in FigOl Decreasing the chemical potential even 
further, the structured dilute liquid persists until /x = — Vi, where there is a phase 
transition between the structured liquid and the gas phase, for < —Vi, the system 
stays empty. 
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Figure 3. Sub- lattices 



2.0.2. If V\ > — 2V2 For > — 2Vi — 2V2 the system is in the dense hquid and for 
/i < —2Vi — 2V2 the stable one is the gas phase. In this case only one phase transition 
at /i = —2Vi — 2V2 is present. 

3. Mean-Field Approximation 

For the mean field analysis, we will follow a sub-lattice strategy j2SI that it is able to 
capture the phase transition that other mean field schemes miss. The chemical potential 
can be rewritten as the sum of four potentials, one for each sub-lattice: 

1 ^ 

^ a=l 

The use of independent chemical potentials may be necessary in a more general 
problem, where each sub-lattice is composed by a different type of molecule. For our 
model we have: 

/il = Ai2 = = /^4, (4) 

which leads, by eq.Q, to 

yu = /i„, a = 1,2, 3, 4. 
Now, rewriting the Hamiltonian , eq.(IT)), in term of the sub-lattices, we get 

'H = -j:Y.H:"iWA)^^^ (5) 

where 

4 

and where the relation between Jij, V\ and V2 is given by: 

Jy2 = J21 = J23 = J32 = J34 = J43 = JiA = J41 = Vi 
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The mean-field approximation here consists in replacing -ff^-^-^ by its mean value, 
given by 

[Kf^) = f,^ + j2T.J^j(^j)^ ^ea. (6) 

13=1 i€l3 

Under this approximation, the density of the sub-lattice /3 is given by p/3 = {aj) and the 
interaction parameter can be writen as 

= H Jij, iea, j e f3, 

what leads to 

V 

(Ha^^) = Pa + Yl ^^i3Pp, i ea, j e (3. 

13=1 

Substituting eq.© in eq.© we get 



(7) 



n 



MF 



- I H ^c^pPp + Poc\<yi + ]^Y.^Yl (^c.f3pl3po 



(8) 

a=l iGa \/3=l 

where the last term corrects for overcounting. Eq.(jHJ is the Hamiltonian in the mean- 
field approximation. With this, we calculate the mean-field grand potential per site, 
that is 

0^^= - keT \n2 



knT 



In cosh 



0=1 



2kBT 



yi ^apPp + Pa 

yf3=l 



(9) 



a=l \/3=l 



a=l 13=1 



where ks is the Boltzmann factor and T is the temperature. 
The density can be derived from this grand potential as: 

/ PiaMF \ 



Po 



what leads to 



Pa 



dpc 



tanh 



a = 1,2,3,4, 



2kBT 



a 



1,2,3,4, 



(10) 



H ea(3p(3 + Pc 

the mean density of a sub-lattice a. 

Solving Eqs. (jTU)) for fixed values of temperature and chemical potential it is possible 
to obtain the density of each sub-lattice and consequently not only the overall density 
of the system but also specific phase. For instance, if for a given temperature and 
chemical potential, the sub-lattice densities would be pi = p4 = and P2 = P3 = 1, the 
system would be in structured dilute liquid phase. Since we know that for T = two 
liquid phases exist only if Vi < — 2V2, we explore this region of the parameter space. 
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Figure 4. Mean field phase diagram : The open circles represent the second order 
phase transition, the solid lines are the two first-order phase transitions and the two 
filled circles are the two tricritical points. 



For simplicity, we set Vi = 1 and V2 = — 1 (other parameter choices will not affect 
qualitatively the result) and we construct the mean-field phase diagram shown in Fig. 



At high temperatures, each sub-lattice is half full in an disorganized way. This is 
the fluid phase. As the temperature is lower at a fixed potential = fi* = fi/Vi > 1, 
all sub-lattices become filled, no phase transition is observed. For very low chemical 
potentials, /i* < —1, as the temperature is decreased, the system goes from the fluid 
to the gas phase continuously with no phase transition. For intermediate chemical 
potentials, —l>fi*>l, there is a continuous phase transition between the fluid phase 
and the structured dilute liquid phase. The coexistence line between the gas phase and 
the structured diluted liquid meets the continuous phase transition at a tricritical point. 
Another similar point is observed at the contact between the coexistence line between 
the structured diluted liquid and the dense liquid and the continuous line. The density 
is monotonic, therefore no density anomaly is observed. 

4. Monte Carlo Simulations 

The rather simple mean field approach introduced in the previous session is unable to 
account for the density anomalies. For investigating the possibility of a density anomaly 
in our potential, Monte Carlo simulations in grand canonical ensemble were performed. 
The Metropolis algorithm was used to study square L x L lattice and IV2I/V1 = 1. 
Different system sizes L = 10, 20, 30, 40, 50 were investigated. The typical equilibration 
time was 1500000 Monte Carlo time steps for each lattice site. 



H 
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Figure 5. Total density :for the 50, 40, 30, 20 and 10 lattices, from top to bottom. 



The nonzero temperature phase-diagram was obtained as follows. For a fixed 
temperature and chemical potential the density of each sub-lattice and the specific 
heat was computed by averaging over 5000 measures. Between consecutive measures, r 
Monte Carlo steps were performed to decorrelated the system. The correlation time r 
was calculated using the density correlation function 

1 tmax t 
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Fig. El shows the total density for the lattice sizes L = 50,40,30,20 and 10, from 
top to bottom (see also inset) and for fixed chemical potential fi* = 0.5. According to 
this graph, one could conclude that no phase transition happens as the temperature is 
decreased at a fixed chemical potential fi* = 0.5. However, analyzing the sub-lattice's 
densities illustrated in Figs. El-El ( here we are illustrating only the result obtained for the 
L = 20), it's clear from the strong density fiuctuations that appear at T* = T/Vi = 0.45 
a phase transition occurs. The peak in the specific heat shown in Fig. ^1 confirms the 
existence of this transition and helps to determined the exact location of the critical 
point T*. 

The energy histograms constructed for fixed temperatures around the critical 
temperature T* are shown in Figs. [TTlfT^ For temperatures above and below T*, 
there is only one peak indicating that the transition is continuous ( two peaks would be 
a signature of first-order transition). 

Analyzing again Figs. El- El it can be seen that for T* > 0.45 the sub- lattice densities 
Pi3 ~ 0.56, indicating that the system is in the fiuid phase. For T* < 0.45, pi = P2 = 
and p3 = p4 = 1, that is a signature of the structured dilute liquid. Hence, it is clear 
that around T* = 0.45 there is a continuous transition between the fiuid and the SDL 
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Figure 10. Specific heat for /j,* = 0.50 and L = 50, 40, 30, 20, 10 from top to bottom. 
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Figure 11. Energy histogram for T < Tc, fi* = 0.50 and L ~ 20 
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Figure 12. Energy histogram for T w Tc, fi* = 0.50 and L = 20 



phases. The error bars associated with the location of the critical temperatures for the 
lattice sizes L = 10, 20 are respectively AT = 4 x 10""^, AT = 7.5 x 10"^ while for the 
lattice size L = 30, 40, 50 is given by AT = 6.5 x 10"^ 

Criticality happens only in the thermodynamic limit. Therefore in order to 
determined the actual critical temperature, it is necessary to extrapolate the results 
obtained for the finite system to the limit of L oo. The critical temperature for 
the finite system can be obtained: from the maximum of the specific heat or from the 
minimum of the fourth order Binder's cumulant Ve |2Z1- This last method requires 
computing for each lattice size, the quantity: 



Ve = 1 

where E is the total energy of the system. 



(12) 
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Figure 13. Energy histogram for T > Tc, n* = 0.50 and L = 20 
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Figure 14. T* x l/i^ for L = 10,20,30,40,50. The circles are the minimum of 
Binder's cumulant and the triangles are the maximum of specific heat. 

Fig. El shows the critical temperature of the finite system for different lattice sizes 
obtained from the two methods refereed above. The extrapolation to 1/L gives 
the critical temperature for the infinite system. The values obtained from the peak 
of the specific heat and from the minimum of the Binder's forth order cumulant are 
respectively 0.467 and 0.466. The difference between the two results is smaller than the 
error on the calculation of both numbers. The lattice L = 10 was excluded from the 
extrapolation because it is so far from the infinite size limit that its inclusion would 
require carrying nonlinear terms into the extrapolation. 

Fig. El shows that for a fixed chemical potential the density has a maximum at a 
certain temperature T^ax- The temperature of maximum density as a function of the 
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Figure 15. /i* x T* phase-diagram : the circles indicates the critical line, the solid 
lines show the first-order transitions, the squares locate the temperature of maximum 
density and the full circles are the tricritical points. The error bars are smaller than 
the symbols. 
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Figure 16. p x /i* for T* = 0.19 



chemical potential is is illustrated on Fig. UHl 

In order to observe the first-order phase transitions between the gas and SDL phase 
and between the SDL phase and the DL phase, simulations at fixed temperature and 
varying the chemical potential were performed. Figure IT?)1 shows that for T* = 0.19 the 
transition between the SDL and the DL happens at /i* = 1. At low chemical potentials 
similar transition is observed between the gas and the SDL phases. At the intersection 
between the critical line and the first-order phase transition lines arise two tricritical 
points as indicated in Fig. [1^1 

By integrating the the Gibbs-Duhem relation: 

SdT -Vdp + Ndfx = 0, (13) 
where S stands for entropy and p for pressure, together with the MC simulations at a 
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Figure 17. p* = p/\Vi\ x T* phase diagram : tlic solid lines indicate first- 
order transitions, the empty circles locate the continuous line, the squares show the 
temperature of maximum density and the filled circles are the tricritical points. The 
error bars are smaller than the symbols. 

fixed temperature, the dependence of tlie pressure on tlie density for a fixed temperature 
was obtained. Comparing tliis dependence for different temperatures, was possible to 
extract tlie pliase-diagram illustrated at FiglTTl Here the density anomaly appears as 
the temperature of maximum density for a fixed pressure. 

The phase-diagram for other values of Vi < — 2V2 is similar to the one we present 
here. In the case of Vi > — 2V2 both the two liquid phases and the density anomaly are 
not present. 

5. Conclusions 

The phase-diagram of a two dimensional lattice gas model with competing interactions 
and in contact with a reservoir of particles and temperature was investigated with the 
purpose of testing this attraction/repulsion potential for the presence of criticality and 
density anomaly . 

This system exhibits two liquid phases and a line of density anomalies. Differently 
from the general belief, the density anomaly line in the present case is not associated 
to a single critical point jHI but with a line of critical points. Besides, the density 
anomaly does not arise from a softened core potential but from an outer shell repulsion 
that competes with a short-range attraction. We believe that the common ingredient 
that leads to the presence of demixing between two liquid phases both in softened core 
potentials and in the present model is the existence of two competing scales for the 
interaction [21! [2H]. 

The relation between the form of the potential, criticality and the density anomaly 
goes as follows. The presence of two liquid phases comes from the presence of competing 
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scales. While the short-range attraction favors the formation of a dense liquid phase, 
the outer shell repulsion favors the formation of an open structure, the SDL phase. 

In systems dominated by short-range attractive forces, if the pressure is kept fixed, 
the density increases on cooling. In our case, similar behavior is only observed at 
high temperatures where short-range interactions are dominant. As the temperature is 
decreased, the outer shell repulsion prevents the density to increase beyond a certain 
limit. Therefore, the same competition responsible for the appearance of two liquid 
phases leads to the density anomaly. 

One should point out that the presence of a critical line instead of a single critical 
point as one could generally expect [0] is not surprising. Due to the lattice structure, 
the SDL is not one single phase but the region where two different phases coexist ( 
empty/full rows and empty/full columns). These two phases become critical together 
with the DL phase at the tricritical point that is also the locus where critical line ends. 
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